Interplay Between Time-Temperature- Transformation and 
the Liquid-Liquid Phase Transition in Water 
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We study the TIP5P water model proposed by Mahoney and Jorgensen, which is closer to real 
water than previously-proposed classical pairwise additive potentials. We simulate the model in 
a wide range of deeply supercooled states and find (i) the existence of a non-monotonic "nose- 
shaped" temperature of maximum density line and a non-reentrant spinodal, (ii) the presence of a 
low temperature phase transition, (Hi) the free evolution of bulk water to ice, and (iv) the time- 
temperature-transformation curves at different densities. 
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Much effort has been invested in exploring the over- 
all phase diagram of water and the connection among its 
liquid, supercooled and glassy states jl], |^, ||, with par- 
ticular interest in understanding the origin of the striking 
anomalies at low temperatures, such as the T-dependence 
of the isothermal compressibility Kt, the constant pres- 
sure specific heat Cp, and the thermal expansivity ap. 

The "stability limit conjecture" attributes the increase 
of the response functions upon supercooling to a contin- 
uous re-tracing spinodal line bounding the superheated, 
supercooled and stretched (negative pressure) metastable 
states [Q. This line at its minimum intersects the tem- 
perature of maximum density (TMD) curve tangentially. 
More recently, a different hypothesis has been developed, 
for which the spinodal does not re-enter into the positive 
pressure region, but rather the anomalies are attributed 
to a critical point below the homogeneous nucleation 
line ||. The TMD line, which is negatively sloped at pos- 
itive pressures, becomes positively sloped at sufficiently 
negative pressures and does not intersect the spinodal. A 
line of first order phase transitions — interpreted as the 
liquid state analog of the line separating low and high 
density amorphous glassy phases || |j| — develops from 
this critical point. 

Simulations of supercooled metastable states are possi- 
ble because the structural relaxation time at the temper- 
atures of interest is several orders of magnitude shorter 
than the crystallization time. It is difficult, but not im- 
possible Q, to observe crystallization in simulations of 
molecular models [Q because homogeneous nucleation 
rarely occurs on the time scales reachable by present day 
computers. Bulk water simulations have been crystal- 
lized by applying a homogeneous electric field || or plac- 
ing liquid water in contact with pre-existing ice || [To| , 
but spontaneous crystallization of deeply supercooled 
model water has not been observed in simulations. 

In contrast, experimental measurements of metastable 
liquid states are strongly affected by homogeneous nucle- 
ation. The nucleation and growth of ice particles from 
aqueous solution has been extensively studied, and the 



"nose-shaped" time-temperaturc-transformation (TTT) 
curves have been measured jl], [H], The non- 

monotonic relation between crystallization rate and su- 
percooling depth results from the competition between 
the thermodynamic driving force for nucleation and the 
kinetics of growth jj] . Crystallization hinders direct ex- 
perimental investigation of pure metastable liquid wa- 
ter below the homogeneous nucleation line; only indirect 
measurements can be made by studying the metastable 
melting lines of ices j| . 

This work attempts to unify the phenomena connected 
with the existence of a liquid-liquid phase transition and 
homogeneous nucleation in a single molecular dynamics 
simulation study. We simulate a system of N = 343 
molecules interacting with the TIP5P potential |l3| |. 
TIP5P is a five-site, rigid, non-polarizable water model, 
not unlike the ST2 model @. The TIP5P potential ac- 
curately reproduces the density anomaly at 1 atm and 
exhibits excellent structural properties when compared 



with experimental data ||13|, |15| . The TMD shows the 
correct pressure dependence, shifting to lower temper- 
atures as pressure is increased. Under ambient condi- 
tions, the diffusion constant is close to the experimental 
value, with reasonable temperature and pressure depen- 
dence away from ambient conditions jl3j . 

We perform equilibration runs at constant T (Berend- 
sen's thermostat), while we perform production runs in 
the microcanonical (NVE) ensemble. After thcrmaliza- 
tion at T = 320 K we set the thermostat temperature 
to the temperature of interest. We let the system evolve 
for a time longer than the structural relaxation time r a , 
defined as the time at which F s (Qo,r a ) = 1/e, where 
-Ps(Qoj t) is the self-intermediate scattering function eval- 
uated at Qo = 18 nm _1 , the location of the first peak of 
the static structure factor. In the time r Q , each molecule 
diffuses on average a distance of the order of the near- 
est neighbor distance. We use the final configuration of 
the equilibration run to start a production run of length 
greater than several t q and then analyze the calculated 
trajectory. We check that no drift in any of the stud- 
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FIG. 1: Dependence on density of the pres- 

sure at all temperatures investigated (T = 
215, 220, 230, 240, 250, 260, 270, 280, 290, 300, 320 K, from 
bottom to top). Each curve has been shifted by n x 150 MPa 
to avoid overlaps. An inflection appears as T is decreased, 
transforming into a "flat" coexistence region at T — 215 K, 
indicating the presence of a liquid-liquid transition. Inset: A 
detailed view of the T = 215 K isotherm. 



ied quantities and no crystallization occurs during the 
production run. 

In Fig. [l] we show results for pressure along isotherms. 
At lower temperatures an inflection develops, which be- 
comes a "flat" isotherm at the lowest temperature, T = 
215 K. The presence of a flat region indicates that a phase 
separation takes place, and we estimate the critical tem- 
perature Tc< — 217 ± 3 K, the critical pressure Pc> = 
340± 20 MPa, and the critical density p c > = 1.13 ± 0.04 
g/cm 3 . 

In Fig. ||(a) we plot the pressure along isochores. The 
curves show minima as a function of temperature; the 
locus of the minima is the TMD line, since {dP/dT)v = 
ap/Kx- It can be seen that the pressure exhibits a min- 
imum if the density passes through a maximum (a = 0) . 
It is clear that, as in the case of ST2 water, TIP5P water 
has a TMD that changes slope from negative to positive 
as P decreases. Notably, the point of crossover between 
the two behaviors is located at ambient pressure, T w 4 
C , and p w 1 g/cm 3 . 

We also plot the spinodal line. We calculate the points 
on the spinodal line fitting the isotherms (for T > 300K) 
of Fig. | to the form P(T, p) = P S (T) + A [p - p s (T)] 2 , 
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FIG. 2: (a) Pressure along seven isochores; the minima cor- 
respond to the temperature of maximum density line (dashed 
line). Note the "nose" of the TMD line at T = 4 C. Stars 
denote the liquid spinodal line, which is not reentrant, and 
terminates at the liquid-gas critical point, (b) The full phase 
diagram of TIP5P water. The liquid-gas critical point C is in- 
dicated by the filled square jl^] and liquid-liquid critical point 
C' by the filled circle. 



where P S {T) and p s (T) denote the pressure and density 
of the spinodal line. This functional form is the mean 
field prediction for P{p) close to a spinodal line. For 
T < 250K, we calculate P S (T) by estimating the location 
of the minimum of P(p). The results in Fig.|| show that 
the liquid spinodal line is not reentrant and does not 
intersect the TMD line @. 

A supercooled liquid is metastable with respect to the 
crystal, so it is driven to crystallize 0]. However, crystal- 
lization of model water has not been found in simulations 
because the homogeneous nucleation time far exceeds the 
CPU time. However, for TIP5P, crystallization times lie 
within a time window accessible to present-day simula- 
tions, and we observe crystallization at densities p = 1.15 
g/cm 3 , and 1.20 g/cm 3 for a wide range of temperatures 
(FigJ). 

To quantify the crystallization process, we analyze four 



3 




10"' 10° 10 1 io 2 io 3 io 4 io 5 

Time [ ps ] 



FIG. 3: Average crystallization time (open circles) as a func- 
tion of temperature at the two densities (a) 1.15 g/cm 3 and 
(b) 1.20 g/cm 3 . A well-defined nose shape is visible, as mea- 
sured for water solutions jl]]]. We also show the structural 
relaxation times r a as calculated from the self-intermediate 
scattering function F 3 (Q,t) (closed circles). Shown as solid 
lines are the MCT power law fits with T c = 211 K, 7 = 2.9 
for p = 1.15 g/cm 3 , and T c = 213 K, 7 = 2.13 for p = 1.20 
g/cm 3 . The interplay between these two curves is discussed 
in the text. 



independent configurations thermalized at temperature 
T = 320 K and instantaneously quenched to the tem- 
perature of interest. We monitor the potential energy 
as well as the time evolution of the structure factor 
S'g(i) = {pQ{t)p*^{t)) /N at all wave vector Q values, 

ranging from the smallest value 2ir/L allowed by the side 
L of the simulation box up to 50 nm _1 . The oxygen den- 
sity fluctuation Pg(t) = Yli=i ex P (*Q ' where is 
the oxygen coordinate of molecule i. The onset of crys- 
tallization coincides with the occurrence of (i) a sudden 
drop in potential energy and (ii) a sharp increase in the 
density fluctuations at one or more wave vector values. 
When crystallization occurs, the value of S^(t) jumps 
from 0(1 ) in the liquid to O(N). 

Defining the crystallization time is somewhat arbitrary 
because of the stochasticity which accompanies the onset 
of crystallization and the definition of the critical nucleus. 
We define T crys t as the time at which any density fluc- 
tuation S^t) grows above a threshold value S* = 15 
and remains continuously above the threshold for a time 




FIG. 4: For p — 1.20 g / cm 3 , the energy minimum configu- 
ration, viewed along the c-axis, of a crystal formed at T = 270 
K. 



exceeding t* = 40 ps. This threshold prevents transient 
density fluctuations from being attributed to crystalliza- 
tion. We also perform calculations for other definitions 
of S* and t*, but the above values are sufficient to un- 
ambiguously identify the onset of crystallization without 
requiring excessive simulation time. Fig. ^ shows the 
crystallization times T cryst , averaged over the four inde- 
pendent runs, for two different densities and for a broad 
range of T. The resulting TTT curve shows a character- 
istic "nose" shape, arising from the competition between 
two effects, the thermodynamic driving force for nucle- 
ation and the kinetics of growth [JIJ |(| . As temperature 
is lowered, both the thermodynamic driving force and 
the relaxation time increase, and it becomes more diffi- 
cult for particles to diffuse to the energetically-preferred 
crystalline configuration. For both densities, p = 1.15, 
and 1.20 g/cm 3 , the T at which nucleation is fastest is 
around 240 K. At this T, the onset of crystallization re- 
quires about 3 ns. At the lowest studied T, the crystal- 
lization time has grown to 30 ns. 
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Fig. H also shows the relaxation times r a . The T de- 
pendence of r a can be described by a power law t q oc 
(T — T c ) 7 , in agreement with the prediction of mode cou- 
pling theory fllTf . Since the relation r Q <C T cryst holds 
at each temperature, including in the deeply supercooled 
region, "equilibrium" studies of metastable water can be 
achieved before nucleation takes place. The liquid can 
be connected to the deeply supercooled state via "equi- 
librium" metastable states if we choose a quench rate 
larger than the critical cooling rate 1Z C = (T m — T n )/r n , 
where T rn is the melting T and T n and r„ locate the 
nose in the TTT curve g. For TIP5P, K c « 10 10 K/s 
at the two studied densities. For p = 1.10 g/cm 3 and 
T = 240 K, we observe only one (out of four) crystalliza- 
tion event within a time of 70 ns. For densities smaller 
than p = 1.10 g/cm 3 , we observe no crystallization events 
within a time of 60 ns and hence we can only estimate 
that 1Z C is smaller than 10 9 K/s (the experimental value 
for water at ambient pressure is 1Z C ss 10 7 K/s flS|). 

In Fig. [I] we show a typical crystal configuration. The 
crystal structure, after energy minimization at constant 
volume, is a proton-ordered structure similar to ice-B, 
first observed by Baez and Clancy pi. Ice-B is a variant 



of the ice IX structure, which is the proton-ordered form 
of ice III. The density of ice IX and ice III is in fact 1.16 
g/cm 3 , close to our value. 

We have shown that the liquid-liquid phase separation 
can be observed in metastable equilibrium (i) if the cool- 
ing rate is faster than 1Z C , and (ii) if the observation 
time is shorter than the crystallization time at the criti- 
cal point. While both such conditions can be realized in 
numerical simulations — as shown here — they cannot 
be met in experiments. Our simulations also show that a 
continuity of states between liquid and glassy phases of 
water exists Q. Liquid states below the crystallization 
temperature can be accessed provided the cooling rate 
exceeds 1Z C . 
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